;procedure that takes an input wavefront and propagates it...
;m_pix: the number of m per pixel
;z:	the propagation distance
;lambda:the wavelength 
;
; Original notes from 29 Oct 02:
;
;I've lost some code that used multiplication by a Fresnel kernel to propagate
;a wavefront. This kernel is the Fourier transform of exp(ik/2z*(x^2+y^2)), 
;which is sqrt(2*!pi/k/z)*exp(-i*u^2*z/2/k), where k is wavenumber, z is the
;propagation distance and u is spatial frequency (in m^-1 or some normalised
;version...). If the code doesn't turn up, then this quick note should enable
;its reproduction.

function fresnel, wf, m_pix, z, lambda

z = double(z)
dimx = (size(wf))[1]
dimy = (size(wf))[2]
if (dimx ne dimy) then begin
 print, 'Input array must be square'
 return, wf
endif

;k = 1.0/lambda;1./lambda ;The wavenumber

ft_wf = fft(wf,1)

;NB Nyquist corresponds to u = 1/2/m_pix
;kernel = sqrt(2*!pi/k/z)*exp(-complex(0,1)*(dist(dimx)/dimx/m_pix)^2*z/2/k)

kernel = exp(-complex(0,1)*(2.0*!pi*dist(dimx)/dimx/m_pix)^2*z*lambda/4.0/!pi)
;kernel = exp(-complex(0,1)*!pi/z/lambda*(dist(dimx)*m_pix)^2)
;kernel = exp(-dcomplex(0,1)*!pi/lambda*( $
;	 (dist(dimx)*m_pix)^2/z )); + (dist(dimx)*m_pix)^4/z^3/4.0 ) )
;kernel = fft(kernel,1)

return, fft(kernel*ft_wf,-1)

end
